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Analytical gradients 1-7 have become a standard tool 
in molecular quantum chemistry. They are indispensable 
for the optimization of structures, and many properties 
can be efficiently computed with the help of analytical 
derivatives. The field was pioneered by Pulay 8 ; the the- 
CIh' ory had already been derived earlier independently 9 . 

The traditional quantum chemical methods are dif- 
' ficult to apply to solids because of the large increase 
, of the computational effort with the system size. Af- 
{SJ ' ter several decades of development, Hartree-Fock calcu- 
lations for solids can nowadays be routinely performed 
with the CRYSTAL code 10 - 11 . Although Hartree-Fock 
calculations often have large errors due to the neglect of 
. electronic correlation, a large interest has grown in the 
past few years due to the success of hybrid functionals 
which include an admixture of exact (Fock) exchange. 
Analytical gradients in the CRYSTAL code were first 
' implemented with respect to nuclear positions 12,13 , and 
after the implementation of a scheme for geometry opti- 
mization, an efficient structural optimization could be 
O |. performed 14 . In periodic systems, the cell parameter 
is another variable to be optimized. The first gradi- 
• i-h , ents with respect to the cell parameter, at the Hartree- 
Fock level, were for systems periodic in one dimension 15 . 
Various groups have implemented these gradients in one 
dimension 16,17 (see also the recent review article 18 ) or 
in two dimensions 19 . For the general case, a strategy to 
compute cell parameter derivatives (and thus the stress 
tensor) was suggested with point charges 20 , and an al- 
gorithm for structural optimization, based on redundant 
internal coordinates was proposed 21 . Second analytical 
derivatives with respect to the cell parameter have also 
been implemented recently 22 . 

A first big step of the corresponding implementa- 
tion in the CRYSTAL code were analytical Hartree- 
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Fock gradients with respect to the cell parameter in 
three dimensions 23 . It is important to note that the 
CRYSTAL code is based on the Ewald 24 ' 25 method in 
three dimensions, so that computing analytical gradi- 
ents with respect to the cell parameter requires vari- 
ous additional derivatives: for example the reciprocal 
lattice vectors depend on the cell parameter, and var- 
ious others. This requires additional derivatives which 
were not yet available with the implementation of nu- 
clear gradients, and this has been documented in great 
detail 23 . The one and two-dimensional case are again dif- 
ferent because different potentials are used: Parry's po- 
tential in two dimensions 26 ' , and Saunders' potential in 
one dimension 28 . Parry's potential is similar to Ewald's 
potential, but modified for the case of two dimensions. 
Saunders' potential relies on a real space approach. 

This article is intended to complement the first article 
on cell gradients 23 . Many parts have already been de- 
scribed in the first article, and therefore the main empha- 
sis is to delineate the differences due to the dimension- 
ality. The article consists thus of one section about the 
general differences to the three-dimensional case, one sec- 
tion about the two-dimensional case, one section about 
the one-dimensional case, and one section with examples. 

II. GENERAL DIFFERENCES WITH RESPECT 
TO THE THREE-DIMENSIONAL CASE 

The main difference to the three-dimensional case is 
the way how the Coulomb energy is computed. The ex- 
pression to be evaluated is the Coulomb energy per cell: 

^coui = 1 J d 3 r J d 3 r'p{r)$(r-r')p{r') 

(1) 

with $ being the potential function corresponding to 
three dimensions (Ewald's potential function) 24 ' 25 , two 
dimensions (Parry's potential function) 26 or one dimen- 
sion (Saunders' potential function) 28 . 

p(r) is a cellular charge distribution, composed of the 
nuclear charges Z a at the positions of the nuclei A a , 

p ™ { 7 )= J2z a S(r-A a ) (2) 

a 

and the electronic charge distribution 

P 6 Vl = - E P ^Mr- A,)Mr- X - g) (3) 

The basis functions cf)^ (r — — g) are real spheri- 
cal Gaussian type functions, -P^g ^ s * ne density matrix 

in real space. A^ denotes the nucleus where the basis 
function fj, is centered. The implementation is done for 
the case of closed shell Hartree-Fock and unrestricted 
Hartree-Fock methods. For the sake of simplicity, the 
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spin is ignored in the equations in this article. The ex- 
tension is straightforward, as was shown for the three- 
dimensional case 23 . Examples for spin-polarized calcula- 
tions are given in section V. 

The potential function enters via the nuclear-nuclear 
repulsion (equation 10 in reference 23), the nuclear at- 
traction integrals (equation 34 in reference 23), and the 
field integrals (equation 43 in reference 23). Essentially, 
the derivatives are computed as described in the previous 
article 23 , there are only minor differences as described in 
section III and IV. 

The derivatives of the other integrals (overlap, kinetic 
energy, multipoles, bielectronics) and the calculation of 
the energy-weighted density matrix is practically identi- 
cal to the three-dimensional case 23 . 

Finally, the correction due to the spheropole (equation 
47 in reference 23) is zero in one and two dimensions and 
thus does not have to be discussed. The spheropole is a 
correction which arises due to the Ewald method, when 
applied to the electronic charge distribution: the charge 
distribution is approximated by multipoles in the long 
range, and not approximated in the short range. The 
electrostatic potential is then computed as the sum of the 
Ewald potential of the multipoles and of the Coulomb 
potential of the charge distribution in the short range. 
Replacing the Ewald potential with the Coulomb poten- 
tial is correct, if the difference of multipolar charge dis- 
tribution and the exact charge distribution in the short 
range, has zero charge, dipolc, quadrupole, and second 
spherical moment 25 . The second spherical moment can 
also be seen as the average electrostatic potential of a 
charge distribution (see the discussion in section 3.2 of 
reference 25). Here, it corresponds to the average elec- 
trostatic potential of the difference of the exact and the 
approximated charge distribution. This term is finite and 
in general non-zero, in the case of periodicity in three di- 
mensions. However, when the system has periodicity in 
less than three dimensions, the average electrostatic po- 
tential of a charge distribution with zero total charge, 
dipole and quadrupol, is zero. Therefore, there is no 
spheropole in less than three dimensions. 

This can also seen from equation 31 in reference 25. 
The average Coulomb potential is obtained as follows: 

* = ~W~ [P diff (rl^d 3 r (4) 

oVcell J 

pdiff ^ corresponds here to the difference between the 
exact charge distribution and the multipolar charge dis- 
tribution. The integral is over the whole space and fi- 
nite. The prefactor involves a division by the cell volume 
V ce ii of the three-dimensional cell. We might now ap- 
proximate a system with periodicity in two dimensions 
by a system of slabs with three-dimensional periodicity, 
where the slabs are separated by a vacuum region. When 
we increase the vacuum region and thus the cell volume 
V ce ii, then the integral remains essentially constant, but 
the prefactor becomes smaller and smaller and therefore 
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the average Coulomb potential becomes zero, and the 
spheropole correction becomes zero. 

It should be mentioned, that two-dimensional period- 
icity is implemented in the CRYSTAL code in such a 
way that there is only one slab which is not repeated in 
the third dimension. Still, the argument presented above 
holds in a similar way, and there is thus no spheropole 
correction in systems with less than three-dimensional 
periodicity 

The total energy is thus similar to the three- 
dimensional case 23 , apart from the spheropole term 
which is zero: 

^total ^kinetic _|_ ^NN _|_ ^coul — nuc _|_ ^coul— cl _|_ ^exch — cl 



E P ^o T »o» S + \ E z aZMA b A a ) 

g.fi.v a.b 

E p ^o E z - J M?- A*)M?- & - g)W- A)d 3 r 

2 E P vgtiO\. E P ahrO^ nOvgrOcrh ~EE E ^ ^P c '^ c ^ M l^ugc) 
g,H,v h.T.a c i=0 m=-l 

o E P vgnO E P ahTQ^- iSvgraah (^) 



9,V,v h,T,<J 

The individual terms contributing to the total en- 
ergy are the kinetic energy E klnctlc , the nuclear- nuclear 
repulsion energy _E NN , the nuclear-electron attraction 
^coui-nuc^ the c i cc tron-electron repulsion E coul -° l an d 
the Fock exchange £: exch - cl . The variables will not all 
be explained in order to reduce the number of formu- 
las in this article. The reader is referred to the article 
on the three-dimensional case for the details where all 
these terms are explained 23 . The gradient with respect 
to the cell parameters a»j is given in the following equa- 
tion. As the total energy, the gradient is similar to the 
three-dimensional case apart from the spheropole term 
which is zero. 

Q Jjtotal 

Fa ' J = 8a tJ = 

"9>» daij datj ^ vS "° datj 



i ( h( - - l. i 

_1 V P JVP- HOvgrOvh \ \^ \ " \ " p^ 

g,H,v T,cr lJ c 1=0 m=-l f i rec a 

^_ [ J Mr- A T )Mr~ A a K)XT{r- A c )d 3 r ] } 

1 r)X 

2 2-~i crhrO Q a . . 
S.M." h,T,a l ° 

dS - [ 

_ / exp(ii?g) V a vn (K)a^ n (K)e n (K)Q(e F - e n (K))d\ (6) 

OQij J BZ 
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III. THE TWO-DIMENSIONAL CASE 



In the two-dimensional case, the primitive cell is given 
by two vectors, with two components: ai, a 2 . a>ij are 
defined in such a way that an = a\ x is the x-component 
of di, ai2 = aiy the y-component of a\, a 2 \ is the x- 
component of a 2 , and a 2 2 is the y-component of a 2 . 



ai \ _ I a lx ai y \ _ I a xl a 12 
0,2 ) \ a 2x a 2y ) \ 021 a 2 2 



(7) 



A point g of the direct lattice is defined as g — n\a\ + 
n 2 a2, with m, n 2 being integer numbers. The position of 
an atom c in a cell at the origin (i.e. g = 0) is given as 
A c = fc.\d\ + / c ,2a2, and then in cell g the position will 
be:^ 

A c + 9 = (fc.i + n SA )a! + (/ Cj2 + n sa )a 2 
We have used an additional index, i.e. means fac- 
tor rii of the lattice vector g. The cartesian t component 
(with t being x or y) of the vector A c + g, indicated as 
A c ,t +9t, is thus 

Act + gt — X)m=l(/ C < m + n g,m)amt 

As all the integrals depend on the position of the nuclei, 
the derivatives of the nuclear coordinates with respect to 
the cell parameters are required: 

— t; = (fc, m + n^ m )S im Sj t = (f Cji + n^i)Sj t (8) 

OCLij 

J m=l 

with the Kronecker symbol Sjt- 

The main difference, compared to the three- 
dimensional case, is Parry's potential function <&(r — A a ) 
that is used: 

Y \r-A a -h\ 



eM2*KK x ( X -A a nK y (y-A a , y ))) _ )} x _ ^ ( , _ , + ^ 

2t/|A-| V V V7 , 




+ exp(-2 7 r|i?|(z-A a , z )) ^1 + erf - - ^ 

-y (* - A a , z )eri(^(z - 4,,*)) |^cxp(- 7 (z - A a , z f) (9) 

ft, are the direct lattice vectors, K the reciprocal lattice 
vectors. V is the area of the two dimensional unit cell, 
7 is a screening parameter which was optimized to be 
7 = (2.4/V 1 / 2 ) 2 , in the two dimensional case. Note that 
this is different to the three-dimensional case 25 where 
7 was chosen as 7 = (2.8/V 1 ^ 3 ) 2 . The prime in the 
direct lattice summation indicates that the summation 
includes all values of the direct lattice vector h, with the 
exception of the case when \r — A a — h\ vanishes. In this 
case, the term — — i — =r- is omitted from the sum. In the 

|i — A a —h\ 
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reciprocal lattice series, the prime indicates that all terms 
with K ^ are included. 

The error function erf is defined as in reference 23, 
equation 12. 

Like the Ewald potential, Parry's potential depends 
on the variables A c , h, V, 7 and K. The derivative 
with respect to the cell parameters thus requires deriva- 
tives with respect to these variables. For the deriva- 
tives with respect to A c and h this is like in the three- 
dimensional case. There are minor changes due to the 
two-dimensionality for the derivatives of the area V, of 
the K- vectors with respect to and of the screening 
parameter 7. 



1. Derivative of the area 

The area is obtained as the magnitude of the cross 
product of the cell parameters: 

V = \Si x a 2 \ = \ai x a 2y - a ly a 2x \ (10) 

If we assume that a\ x a 2y — a\ y a 2x is positive, then the 
derivatives are obtained as: 



dV 


(11) 


— a 2y 
dai x 


dV 


(12) 


oai y 


dV 


(13) 


— a iy 
oa 2x 


dV 


(14) 


Q = a\ x 
oa 2y 



Essentially, the formulas for the three-dimensional case 
can be used, when setting a\ z = 0, a 2z = and 
^3 = (0,0, 1). This holds also for the derivatives of the 
reciprocal lattice vectors, as described in the following 
paragraph. 



2. Derivative of the reciprocal lattice vectors 
The reciprocal lattice vectors K can be expressed as 
/f = n 1 fe 1 +n 2 62 (15) 

with the primitive vectors bi of the reciprocal lattice 
defined as: 

01 = y (a 2y ,-a 2x ) ; b 2 = — {-a ly , a lx ) (16) 

The derivatives are thus: 

<%i dV 61 



da lx da lx V 



(17) 
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and 



dV h 



da\y daiy V 



(18) 



da 2 



da 2x V V 



(19) 



dh 



■2y 



dV b 1 2tt . . 

■^—^ + ^1,0 

da 2y V V 



(20) 



d6 2 



dV b 2 2tt ,„ „ , 

-75—17 + 77 M 



(21) 



^ W 62 + £(-1,0) 



daiy 



daiy V V 



(22) 



db 2 



dV b 2 



da 2x da 2x V 



(23) 



db 2 



da 



2y 



dV b 2 
da 2y V 



(24) 



3. Derivative of the screening parameter 

The derivative is straightforward, like in the three- 
dimensional case: 

dj_ = dj_ dV_ = _^dV_ 

daij dV da>ij V da^ 

As a whole, Parry's potential leeds to similar terms 
appearing in the derivatives as in the case of the Ewald 
potential. This is what was to be expected, as Parry's 
potential is essentially obtained when Ewald's approach 
to treat the Coulomb interaction is applied to a system 
with two-dimensional periodicity. 



IV. THE ONE-DIMENSIONAL CASE 



In the one-dimensional case, there is only one cell pa- 
rameter: a xx — an = a = \a\. This case is somewhat dif- 
ferent from the two- and three-dimensional case because 
a pure real space approach is used in the CRYSTAL code 
for the potential to describe the Coulomb interaction 28 . 
The potential consists of a point charge +1, neutralized 
by a uniform charge distribution of length a, with charge 
density — K The uniform charge distribution is then 
again compensated. Up to a certain range, the sum- 
mation is performed exactly. For larger distances, the 
summation is instead approximated with the help of the 
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Euler-MacLaurin summation rule. As a whole, the fol- 
lowing expression was obtained 28 : 

. 1 H(U — z,a) + H(U + z,a) _. 

*(^= £ ' |f~^gT o ^+£(M,r)+£(M,-r) (26) 

n=-M ' ' 

The first term comprises the exact part, the next two 
(with the H function) the region due to the uniform 
charge density in the range of the exact sum (from —Ma 
to Ma), the remaining two terms (the ^-function) are the 
approximated part. The prime indicates that terms with 
\r — na\ = are omitted. M is thus the number of cells, 
where the sum is performed exactly, and U — a(M + i). 
a is defined as a = x 2 + y 2 , with r = (x,y,z). H is 
the function H(p,a) = ln(\/p 2 + a + p). £(M,r) and 
£(M, —r) are contributions from the long range part, 
which is approximated by the Euler-MacLaurin rectange 
rule summation formula. For more details, see reference 
28. For the present purpose, it is important to note that 
the direct lattice vector a appears in the potential, but 
no screening parameter 7 and no reciprocal lattice vec- 
tors K as in the two- and three-dimensional case. This 
means that derivatives with respect to the nuclear co- 
ordinates A c and derivatives with respect to the direct 
lattice vectors na appear, which are essentially given by 
the nuclear gradients, multiplied with the fractional co- 
ordinates. The derivatives with respect to a due to the 
H and £ function are very lengthy, but still straightfor- 
ward. They are thus not discussed here, but formulas 
can be derived from Saunders' article 28 . 



V. EXAMPLES 

In this section, we give some numerical examples of the 
accuracy of the gradients. The tests considered are es- 
sentially identical or similar to the test cases distributed 
with the CRYSTAL code and with the ones from ref- 
erence 14. Note that the fractional coordinates of the 
atoms were not optimized. 

First, two systems with one-dimensional periodicity 
are considered. In table I, SN is periodically arranged. 
The analytical and numerical derivative agree well up to 4 
digits, and the minimum of the energy at a=4.42 A agrees 
with the place where the gradient changes its sign. In ta- 
ble II, such a comparison is done for polyglycine. The 
agreement of numerical and analytical gradients is simi- 
lar to SN, and again the vanishing of the gradient agrees 
with the minimum of the energy, to at least 0.01 A. In 
table III, ferromagnetic NiO is studied at the level of 
unrestricted Hartree-Fock. The agreement of numerical 
and analytical gradient can be improved by increasing the 
'TTOL"-parameters 11 , as described earlier 12 ' 23 . Indeed, 
when increasing them from default values to higher ones, 
symmetric in ITOL4 and ITOL5, then analytical and nu- 
merical gradient match better. Note that, when running 
at lower ITOL parameters, an inaccuracy is introduced 
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in the total energy expression and thus in the numerical 
gradients as well. The fact that numerical and analyti- 
cal gradients match less well at low ITOL values is thus 
a combination of an inaccuracy in the energy expression 
(which affects the numerical gradient) and an inaccuracy 
in the analytical gradient. Still, in all the tests performed 
so far, no severe error was found when using default val- 
ues for the ITOL parameters. Using higher ITOL pa- 
rameters is mainly useful for tests of the correctness of 
the code. 

Then, various systems with two-dimensional periodic- 
ity are considered. In tabic IV, 3 MgO layers are con- 
sidered. Numerical and analytical derivative agree to 3 
digits, and the minimum of the energy and the vanish- 
ing of the gradient agree also well. The same accuracy is 
found for AI2O3 in table V, where a slab with 6 atomic 
layers is considered. In table VI, a &2O3 slab was chosen 
as an example for unrestricted Hartree-Fock. The accu- 
racy is slightly worse when comparing the numerical and 
the analytical gradient. This can again be improved by 
increasing the "ITOL" -parameters. The minimum in the 
energy agrees already with default "ITOL" values to at 
least 0.01 A. Finally, in table VII, LiF was arranged with 
two dimensional periodicity, without symmetry, in such 
a way that three components of the cell gradient (a\ x , 
aiy, a,2y) can be computed independently. This test thus 
demonstrates that these components are correctly com- 
puted. 

In table VIII, the CPU times are displayed. The cal- 
culations were performed on a single CPU of a Compaq 
ES45, with a clock rate of 1 GHz. As in the three- 
dimensional case, we compare again the CPU time for 
the integrals with the time for the gradients. The CPU 
time for all the gradients (nuclear and cell gradients) is 
roughly five to ten times the CPU time for the integrals. 
This may become smaller in the future with further opti- 
mizations in the gradient code. Note that the CPU time 
for the self consistent field calculations is relatively high 
because a very low convergence threshold was chosen in 
order to ensure the accuracy of the succeeding gradient 
calculation (the gradient calculation is the more accurate, 
the more accurately the self consistent field equations are 
solved) . 

The CPU times thus indicate that analytical gradients 
can be computed at a relatively low expense. Compared 
with numerical gradients, it appears that analytical gra- 
dients should usually be favorable, especially because nu- 
merical gradients will depend on the step size, and often 
it will be necessary to break a symmetry for a finite dis- 
placement, to compute the numerical gradient. Numeri- 
cal gradients require at least one additional energy eval- 
uation for each coordinate to be optimized, which makes 
analytical gradients clearly favorable, if there is a large 
number of geometrical parameters. 
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VI. CONCLUSION 



A formalism for the calculation of the analytical gra- 
dient of the Hartree-Fock energy, with respect to the cell 
parameter, has been presented and implemented in the 
code CRYSTAL, for the case of systems periodic in one 
and two dimensions. The implementation includes the 
cases of spin-restricted and unrestricted polarization. It 
was shown that a high accuracy can be achieved. 
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TABLE I. SN, with one-dimensional periodicity. A comparison of analytical and numerical gradient is done for various unit 
cell lengths. A [3s2pld] basis set was used for S, and a [2s lp] basis set for N. 



a analytical derivative numerical derivative energy 

[A] [E h /a ] [E h /a Q ] [E h ] 

4.30 0.04144 0.0414 -893.870081 

4.41 0.00372 0.0037 -893.874639 

4.42 0.00064 0.0006 -893.874680 

4.43 -0.00238 -0.0024 -893.874663 
4.500 -0.02208 -0.0221 -893.873013 



TABLE II. Polyglycine. A comparison of analytical and numerical gradient is done for various unit cell lengths. Basis sets 
of the size [2slp] were used for C, O, N and a [Is] basis set for H. 



a analytical derivative numerical derivative energy 

[A] [E h /a ] {E h /a Q } [E h ] 

7.30 0.01956 0.0196 -408.220173 

7.42 0.00116 0.0012 -408.222495 

7.43 -0.00030 -0.0003 -408.222503 

7.44 -0.00175 -0.0017 -408.222484 
7.50 -0.01018 -0.0102 -408.221807 



TABLE III. NiO, ferromagnetic, unrestricted Hartree-Fock. The gradient with respect to the cell parameter is computed for 
two different values of the ITOL parameters. A [5s4p2d] basis set for Ni was used, and a [4s3p] basis set for O. 

a analytical derivative numerical derivative energy 

[A] [E h /ao] [E h /ao] [E h ] 

ITOL 6 6 6 6 12 (default) 

5.00 -0.10864 -0.1074 -1581.454974 

ITOL 6 6 6 12 12 

5.00 -0.10782 -0.1078 -1581.456358 



TABLE IV. MgO surface, 3 atomic layers. The unit cell consists of 3 Mg and 3 O atoms, with ai x = a-2 y = a. Basis sets 
: the size [ 
contribute. 



of the size [3s2p] were used. The derivative with respect to — ® dais _|_ *9 ^ a ^v j g displayed, Jt an d ^ do not 

da oai x da oa^y da da\ y da2 X 



a analytical derivative numerical derivative energy 

[A] [E h /ao] [E h /a ] [E h ] 

2.80 0.10544 0.1058 -823.930493 

2.88 0.01035 0.0108 -823.939034 

2.89 0.00006 0.0006 -823.939142 

2.90 -0.00991 -0.0095 -823.939058 
3.00 -0.09403 -0.0937 -823.928906 
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TABLE V. AI2O3, 6 atomic layers. The unit cell consists of 6 Al and 4 O atoms, with ai x = y/Z/2 * a 2y = \/3/2 * a. Basis 

sets of the size [3s2pldl for Al and [2slpl for O were chosen. The derivative with respect to = ^ dai x ^_ d da2y . g 

da aai x da da 2y da 

displayed, — and — do not contribute. 

oaiy da 2x 

a analytical derivative numerical derivative energy 

[A] [E h /ao] [E h /ao] [E h ] 

4.20 0.27548 0.2757 -1400.244182 

4.40 0.00590 0.0059 -1400.295000 

4.41 -0.00570 -0.0060 -1400.295003 

4.42 -0.01712 -0.0171 -1400.294787 
4.70 -0.27847 -0.2786 -1400.211859 



TABLE VI. Cr203, 6 atomic layers, ferromagnetic, unrestricted Hartree-Fock. The unit cell consists of 6 Cr and 4 O atoms, 
with ai x — V3/2 * a 2y = x/3/2 * a. Basis sets of the size [5s4p2d] for Cr and [3s2p] for O were chosen. The derivative with 



d 

respect to — 
oa 



d dai x d dai y 



dai x da da 2y da 



is displayed, 



d 



da 



and 



d 
da-z 



do not contribute. 



a 


analytical derivative 


numerical derivative 


energy 


[A] 


[E h /ao] 


[Eh/ao] 


[Eh] 






ITOL 6 6 6 6 12 (default) 




4.70 


0.13465 


0.1379 


-4622.589785 


4.87 


0.00426 


0.0069 


-4622.612278 


4.88 


-0.00253 


0.0001 


-4622.612339 


4.89 


-0.00921 


-0.0066 


-4622.612277 


5.00 


-0.07676 


-0.0745 


-4622.603638 






ITOL 6 6 6 12 12 




4.88 


-0.00116 


-0.0011 


-4622.617935 


5.00 


-0.07539 


-0.0754 


-4622.609006 



TABLE VII. LiF, with a unit cell of ai x = 5 A, a 2y = 4 A, and an angle of 60°, resulting in ai y = 2.5 A. The F atoms are 
at (x=0.1, y=0 (x and y in fractional units), z=0.1 A), (x=0.5, y=0.5 (x and y in fractional units), z=0.3 A), the Li atoms at 
(x=0.5, y=0 (x and y in fractional units), z=0.2 A), and (x=0, y=0.5 (x and y in fractional units), z=0.4 A). A [2s lp] basis 
set was used for Li, a [4s3p] basis set for F. 

component analytical derivative numerical derivative 

lE h /a ] [Eh/ao] 



dE 
daix 

dE 
da\ y 

dE 
da 2y 



0.04045 0.0406 



-0.04415 -0.0441 



-0.01838 -0.0183 
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TABLE VIII. CPU times for one single point calculation of the various systems. The calculations were performed on a 
Compaq ES45, using a single CPU (1 GHz). The CPU times refer to the part for the integrals (all the integrals were written 
to disk), the self-consistent field (SCF) procedure, and to the calculation of all the gradients (i.e. nuclear gradients and cell 



gradients). 


system 




CPU time, in seconds 






integrals 


SCF 


gradients 


SN 


1 


1 


6 


Polyglycine 


2 


4 


17 


NiO 


2 


14 


9 


MgO 


5 


3 


52 


AI2O3 


8 


12 


78 


Cr 2 3 


27 


153 


176 


LiF 


3 


18 


20 
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